home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / pskysl.z / pskysl
Text File  |  1996-03-14  |  10KB  |  331 lines
  1.  
  2.  
  3.  
  4. PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))                                                          PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DSKYJ, FBSJ, DSKYRJ, FBSRJ - parallel skyline symmetric linear system
  10.      solver
  11.  
  12. DDDDEEEESSSSCCCCRRRRIIIIPPPPTTTTIIIIOOOONNNN
  13.      PSKYSL solves skyline, or profile, symmetric linear systems of the form
  14.         Ax = b
  15.      where A is an n x n symmetric input matrix stored in either the Jenning's
  16.      profile format or the reverse Jenning's profile format (described later),
  17.      b is an input vector of length n, and x is an unknown vector of length n.
  18.      PSKYSL uses a direct method: A is factored into the form
  19.        A = L D L-transpose
  20.      where L is a lower triangular matrix with unit diagonal and D is a
  21.      diagonal matrix.
  22.  
  23.      The PSKYSL library contains four routines.  DSKYJ() factors the matrix A,
  24.      which is stored in the Jenning's format, into L and D.  FBSJ() solves for
  25.      a vector x using the results (L and D) that are generated by DSKYJ(),
  26.      given an input vector b.  DSKYRJ() factors the matrix A, which is stored
  27.      in the reverse Jenning's format, into L and D. Similar to FBSJ(), FBSRJ()
  28.      solves for a vector x in the reverse Jenning's case, for the given input
  29.      vector b.  It is assumed that the input matrix A has been reordered by a
  30.      popular heuristic such as Gibbs-King, or Reverse Cutill-McGee, for
  31.      profile size reduction before a factorization routine is called.  The
  32.      PSKYSL solvers also have the ability to do incremental factorization of
  33.      the input matrix A so as to reduce the memory requirements to run the
  34.      solvers.
  35.  
  36.      To illustrate the skyline storage schemes, consider the following
  37.      symmetric 4x4 array:
  38.  
  39.              1.0 0.0 2.0 0.0
  40.                  3.0 0.0 4.0
  41.                      5.0 0.0
  42.              symmetric   6.0
  43.  
  44.      The Jenning's representation of the array in FORTRAN would be:
  45.  
  46.              pointers[] = {1, 2, 5, 8}
  47.              values[] = {1.0, 3.0, 2.0, 0.0, 5.0, 4.0, 0.0, 6.0}
  48.  
  49.      And the reverse Jenning's representation would be:
  50.  
  51.              pointers[] = {1, 2, 3, 6, 9}
  52.              values[] = {1.0, 3.0, 5.0, 0.0, 2.0, 6.0, 0.0, 4.0}
  53.  
  54.      Due to symmetry, only the upper triangular part of the array needs to be
  55.      stored.  Each column is stored in a top-down/bottom-up fashion, i.e.
  56.      from the first nonzero entry to the diagonal or the other way around,
  57.      into values[].  The zero entries outside the skyline profile are not
  58.      stored.  The book-keeping array pointers[] keeps track of the locations
  59.      of the diagonal elements in values[].  Note that in the reverse Jenning's
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))                                                          PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      case, the last entry in pointers[] provides the information of the total
  75.      number of terms plus one.
  76.  
  77.      The number of processors that are used for the numerical factorization
  78.      are controlled through an argument that is passed to the factor routines.
  79.      No special environment variable setting is required for this matter.
  80.      Parallel threads (if any) in the PSKYSL library are created with 'sproc'
  81.      system calls, and are not inherited from existing C or Fortran threads.
  82.      Therefore, it is advisable to terminate any outstanding Fortran or C
  83.      slave threads with mp_destroy, etc., before the PSKYSL library is called.
  84.  
  85.  
  86. FFFFOOOORRRRTTTTRRRRAAAANNNN SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  87.      SUBROUTINE DSKYJ (VALUES, POINTERS, DINV, NCPU, NEQ, LEAD)
  88.  
  89.          DOUBLE PRECISION VALUES( * )
  90.  
  91.          INTEGER POINTERS( * )
  92.  
  93.          DOUBLE PRECISION DINV( * )
  94.  
  95.          INTEGER NCPU
  96.  
  97.          INTEGER NEQ
  98.  
  99.          INTEGER LEAD
  100.  
  101.      SUBROUTINE FBSJ (VALUES, DINV, N, POINTERS, B, SCR)
  102.  
  103.          DOUBLE PRECISION VALUES( * )
  104.  
  105.          DOUBLE PRECISION DINV( * )
  106.  
  107.          INTEGER N
  108.  
  109.          INTEGER POINTERS( * )
  110.  
  111.          DOUBLE PRECISION B( * )
  112.  
  113.          INTEGER SCR ( * )
  114.  
  115.  
  116.  
  117.      SUBROUTINE DSKYRJ (VALUES, POINTERS, DINV, NCPU, NEQ, LEAD)
  118.  
  119.          DOUBLE PRECISION VALUES( * )
  120.  
  121.          INTEGER POINTERS( * )
  122.  
  123.          DOUBLE PRECISION DINV( * )
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))                                                          PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))
  137.  
  138.  
  139.  
  140.          INTEGER NCPU
  141.  
  142.          INTEGER NEQ
  143.  
  144.          INTEGER LEAD
  145.  
  146.      SUBROUTINE FBSRJ (VALUES, DINV, N, POINTERS, B, SCR)
  147.  
  148.          DOUBLE PRECISION VALUES( * )
  149.  
  150.          DOUBLE PRECISION DINV( * )
  151.  
  152.          INTEGER N
  153.  
  154.          INTEGER POINTERS( * )
  155.  
  156.          DOUBLE PRECISION B( * )
  157.  
  158.          INTEGER SCR ( * )
  159.  
  160.  
  161.  
  162.  
  163. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  164.      values (input to DSKYJ and DSKYRJ)
  165.              The upper triangular part of the input matrix A stored in either
  166.              the Jenning's or the reverse Jenning's format.
  167.  
  168.      values (output from DSKYJ and DSKYRJ)
  169.              The factor L-transpose, except that the diagonal locations are
  170.              replaced with D, stored in the Jenning's or the reverse Jenning's
  171.              format.
  172.  
  173.      values (input to FBSJ and FBSRJ)
  174.              The factor L-transpose, except that the diagonal locations are
  175.              replaced with D, stored in the Jenning's or the reverse Jenning's
  176.              format.
  177.  
  178.      pointers (input)
  179.              An integer vector pointing to the position of the diagonal
  180.              elements of the input matrix A in values[].  This vector is
  181.              usually of length n in the Jenning's case, and n+1 in the reverse
  182.              Jenning's case.  However, in an incremental factorization, the
  183.              length may be less than n/n+1.
  184.  
  185.      dinv (input to DSKYJ and DSKYRJ)
  186.              A double precision vector of length n.
  187.  
  188.      dinv (output from DSKYJ and DSKYRJ)
  189.              The inverse of D.
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))                                                          PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      dinv (input to FBSJ and FBSRJ)
  207.              The inverse of D.
  208.  
  209.      ncpu (input)
  210.              The number of processors that are used for the numerical
  211.              factorization routines DSKYJ and DSKYRJ.
  212.  
  213.      neq (input)
  214.              The number of equations to be factorized in the call.  The neq
  215.              number is set to n if no incremental factorization is performed.
  216.  
  217.      lead (input)
  218.              The first equation to be factorized in the call.  The lead number
  219.              is set to 1 if no incremental factorization is performed.
  220.  
  221.      n (input)
  222.              The number of rows and columns in the matrix A, or the order of
  223.              A.  n >= 0.
  224.  
  225.      b (input)
  226.              The right-hand-side vector in a FBSJ/FBSRJ call.
  227.  
  228.      b (output)
  229.              The solution vector in a FBSJ/FBSRJ call.
  230.  
  231.      scr (input)
  232.              An integer scratch vector of n long.
  233.  
  234.  
  235. TTTTUUUUNNNNIIIINNNNGGGG
  236.      Optimized and parallelized for the SGI R8000 platform.
  237.  
  238.  
  239. IIIINNNNCCCCRRRREEEEMMMMEEEENNNNTTTTAAAALLLL FFFFAAAACCCCTTTTOOOORRRRIIIIZZZZAAAATTTTIIIIOOOONNNN
  240.      Sometimes it is desired to perform the factorization of the input matrix
  241.      A in an incremental fashion so that a large portion of the array can be
  242.      kept in a disk device to conserve memory resources.  The PSKYSL library
  243.      provides such a memory-saving capability if a proper sequence of
  244.      values[], pointers[], and other book-keepings is constructed and passed
  245.      to a series of calls to the factorization routines.  The following
  246.      example illustrates a single thread incremental factorization of the a.m.
  247.      4x4 array with two calls to DSKYJ.
  248.  
  249.      1) First call
  250.  
  251.              pointers0[] = {1, 2, 5}
  252.              values0[] = {1.0, 3.0, 2.0, 0.0, 5.0}
  253.              dinv0[] = { 0.0, 0.0, 0.0}
  254.              ncpu0 = 1
  255.              neq0 = 3
  256.              lead0 = 1
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))                                                          PPPPSSSSKKKKYYYYSSSSLLLL((((3333FFFF))))
  269.  
  270.  
  271.  
  272.      2) Second call
  273.  
  274.              pointers1[] = {1, 4, 7}
  275.              values1[] = {values0[3], values0[4], values0[5],
  276.                          values0[6], 4.0, 0.0, 6.0}
  277.              dinv1[] = { dinv0[2], dinv0[3], 0.0}
  278.              ncpu1 = 1
  279.              neq1 = 3
  280.              lead1 = 3
  281.  
  282.      Note that parts of the 'values' and the 'dinv' arrays in the second call
  283.      inherit the results from the first call.
  284.  
  285.  
  286.  
  287.  
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.